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I'HK fJAm’ ANALYTIC ALiiODYNAMIC OiLVG .MKTIIOD 


FOK Tffi: VlNTf SATCLLITF TIIKOHV 

S. Watson 
G. D. Mistretta 
N. L, Bonavito 


ABSTilACT 


In order to retain separability in the Vinti theory of earth satellite motion 
when a non-consexwative force such as air dra^ is consi<lercd, a set of variational 
equations for the orbital elements are introduced, aiid expressed as functions of 
the transverse, x'adial, and normal components of the non-conservative forces 
acting on the system. In this approach, the Hamiltonian is preser\^ed in form, and 
remains the total energ>', but the initial or boundary conditions and hence the 
Jacobi constants of the motion advance with time through the variational equations. 
In particular, the atmospheric density profile is written as a ’fitted’ ex]xonential 
function of the eccentric anomaly, which reproduces tabular values of static model 
atmospheric densities at all altitudes to within ninety-eight percent and simulta- 
neously reduces the variational equations to indefinite integrals with closed 
form evaluations, whose limits are in terms of the eccentric anomaly. The values 
of the limits for any arbitrary time interval are obtained from the Vinti program. 

Results of the BMW (Bonavito, Mistretta, Watson) theory for the case of the 
Intense air drag satellites San Marco-2 and Air Force Cannonball are given* 
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These results indicate that the satellite ephemerldes produced by the 
th'fOiy in conjunction with the Vinti program are of very liigh accuracy. In ad~ 
ditiony since the program is entirely analytic, several months of cphemerides 
can be obtained within a few seconds of computer time-. 
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THK liMXV analytic: AKUODYNAMK: DiUC'. MKTIIOl) 
rou Tnk \ inti sathlu i l tiikoky 

IN'I'RODLJCTION 

Jn this paper we shall treat the variation of Izsak eleinents caused by dis- 
turbing forces actin;^ on an orbit, 1'hese elements, an intrinsic part of the \'inti 
theoiy, are obtained as functions of time and then by a process of re-initialization, 
the Vinti equations of motion are solved for the position and velocity of the satel- 
lite. We should note tliat the air resistance of a body has in |j;eneral six com- 
ponents, v.liL'co being forces and three being moments of forces, which tend to 
make the motion of an asymmetric body very complex. However, for the case of 
a non-rotating sphere, the resistance cim be reduced to a single component di- 
rected oppositely to the velocity of the sphere. In this treatment we shall intro- 
duce senii-empiricaliy determined components of drag force which are functions 
of the eccentric anomaly. With use of these, the differential equations for the 
variation with time of che Izsak elements can be solved along with the Vinti equa- 
tions of motion. 

I. DRAG variational EQUATIONS 

Here we present equations for the time rate of change of the Izsak orbital 
elements of the Vinti satellite theory (Reference 1). h\ our treatment, we shall assume 


1 


that these equations represent instajilaneous departures due to aercKiynamic dr:q; 
or perturbinj 5 forces on a Vinti orbit. These depai*tures in the Vinti-Izsak or!)ital 
elements ai'e then integrated as functions of eccentric :uiomaly for the InsUuit of 
time under consideration. 

The differential equations representing valuations of orbital elements ai'c: 
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where n is related to a by - GM^ == //. G Is the gravitational constant, and 

Mg is the mass of the Earth. Any perturbing force F can be represented 



instanlanoously as a potential Ki*adicnt. I.ot i; l>o that pote ntial. We now 
resolve the perturbing; force F into the following; components: 


R, in the direction of the position vector from the force center to the satellite. 

T, is the perpendicular to R, lies in the orbital pkme, and is positive in the 
direction of motion. 

W, is mutually perpendicular to U and T and completes a rignt handed set of 
component directions* 


The partial deruatives appearing in equations (1) are given by: 
.'^U fv\ 

R (it) 

(^V ^ / 2 * o cos V \ 

^ - K a cos V r T I — j — j r sin V 


i!}L 

rV«T 


r W s i ti 0 


rVU 


T r cos i - W( r cos sin i ) 




r T 






( 2 ) 


3 



whex'c is the argunjcnt of latitude. Suhstituting (2) h)to (1), 
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The Vinti parameter related to the orbit inclination is 
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From this, we have 
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or, 


• IS 


2 /S / 1 - S 
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The variation of the inclination can then !)o vvrittcn a« 
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The mean anomaly is related to the time of perij^ec passage ^ by the expression, 
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together with dM/dt ^‘rom (3), we find 



Summarizing our equations for the time rate of change of orbital elements, wo 
have: 
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Let us now substitute for R, T, and W in terms of the eccentric anomaly as 
given in refei'ence 2, page 165: 
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where V is the velocit^^ of the satellite doiined by the equation 
V (1 ^ c cos K )“ ^ ^ ( I - <-• c:os E ^ “ j^( 1 ^ t os E ) 


- «U 1 > (• 
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d . ^ : r '/2 . ,. 2 ) 1/2 (1 . S ) l /2 

where 

Ci)^ = angular velocity of rotation of earth 
fj = atmospheric density 
A -- projected area of satellite 
m = mass of satellite 
Cjj = drag parameter (usually around 2.2) 
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sin V 


(1(5) 
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Since the coordinate 77 " (P ^ Q sin /') sin i , and the quanUty * of tlie Vinti 
potCMitial is approximately 7 kilometers, we then assitmc r ' ‘ 1 - c cos E), 

where .> is an oblate spheroidal coordinate,. Now substituting equations (14) and 
(15) into equations (13), and inserting these residts with equations ( 10 ) into 
equations ( 12 ), we get 
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If we indicate the eccentric anomaly at time t ^ by Kj and at time t^ by , the 
solutions of equations (17) are given by the following: 
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U. ATMaSPIIKRIC DKNSITY KPPUKSPNTATION 

At this point we consider the expression for atmospheric density variation 
given by King-IIele (Reference 3). An expression for the atmospheric density 
as a function of the eccentric anomaly is given by 


p = {1 + b (r - i f-'xp 


= {1 + bx^ (1 - cos E)^} 



cos 



( 1 ) 


where x = ae, = density scale height at perigee which in equation (1) is 
assumed to vary linearly with the perigee height. 

Let us r Av write equation (1) in a more general form so as to represent 
the atmospheric density variation over an^' interval within the orbit in which 
the boundary conditions are known. A generalization of (1) is desirable for 
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our purposes to insure that tho resultant analytic atmospheric moclcl rif^iclly 
adheres to a tabular set of densities at all altitudes. To obtain this generalization, 
let the subscripts L and U refer to a lower and upper point of an interval 
respectively. This interval is judiciously chosen (by a method discussed 

later in this section) to allow /< to be given by an expression of the type of 
equation (1), over this restricted domain [E^^, E^,l, Let, 


P(r)=pL I' ~ • ‘■j, - - ‘‘u 

= {1 + b(L, U) x^(c'os fosE)^}pxp ^-^(cosEj -cosE)|, £| <E<E„ 


( 1 . 1 ) 

where available from density tables, and are available 

from the orbit theory (as are r^ and r^), and b(L, U) is derivable from (1.1) 
by setting p = p^ and forming the inverse transformation 

l>(L. U) r. b(L, U; p^. Py. r^. H^). (1-2) 

If equation (1.1) above is substituted into (1-18), the resulting expressions 
for the variational equations reduce to the form 

j sin" E cos*” E exp {- (x/H^) (cos E^ - cos E)} dE (2) 
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where is held edrstant over Uie ranv;c ol a suhinlcrval 61 integration and x is 
held constant over the entire range of integratioru Since however, the ijaramoter 
n can assume the value zero, we ai'e forced to alter the form of (2) so as to insure 
that the expression under the integral sigh is integrable while simultaneously 
maintaining its theoretical physical consent. In short, the variational equations 
cannot he evaluated by the Fundamental Theoz^em oi' Calculus since an antideri- 
vative of sin" E cos"‘ E exi^ ' - (x/lJj^) (cos - cos E): is not expressible in 
terms of elementary functions. 

To achieve this purpose, let us now consider the following; 

Rewrite the expression for atmospheric density (1) in the form 

[1 + b(L. U) x^(cos - cos E)^l b, (E) exp [b^fEVE. 

< E < E^ (3) 

where bj (E)> b 2 (E) are arbitrary' functions of E which w’ill be chosen to preserve 
closed form integrability of the variational equations (1-18) while retaining near 
complete numerical agreement with the exponential term of the King-Ilelehs 
expression for (E). 

The selection of the two functions bj(E), b^(E) between an upper limit of 
integration ^ lower limit Ej will be made in the following manner, . Con- 

sider the Interval [E^ , E*] subdivided into a collection of n nonoverlapping sub- 
intervals I , = 1, 2, . , • , n), defined by the paxtitioning 
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For any subinterval I = [x. t x ), tlefine 

"■ - 1 


b, (E) = c for X. <E<x. 
' *: (.-I 


H) 


bj (E) = d 


where c^, d, are constants. The determination of c^, dj will be accomplished 
by applying the technique of classical weighted least squares. In this manner 
the definition of bj (E) and b^ (E) over the domain [E*, E*1 are the step 
functions 


bj(E) = Cj 


^ E < Xj 
< E < Xj 


= c X - ^ E < X 

n n 1 n 


bj(E) = d, < E < X, 


= dj X, 2 E < Xj 


( 5 ) 


= d X , < E < X . 

n n 
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Consider observational data Y- j^enerated with 1/ by repcjatec! evaluation of 
the function 

li(Ej) ~ oxp [-'i(cos E, cos E^)i 

whore G = x/ll^ and < E. < Xj. To define the pseudo- rcj;ression situation 
governing our estimation problem, consider the implicitly linear, exponential model 

Y. = /3j €. exp (/^2 ^ i n I, 2, • • • , m (0) 

relating the concomitant variable E and the dependent or z'esponse variable Y 
that, luilike the general regression response variable, displays no random variability 
Hence, the random variables t . denoting the dispersion cnaracteristics of Y% 

about the theoretical regression lino become meaningless in the curve fitting 
analysis. That is, the random varLablos . are degenerate in the sense that 
their probability density functions p ) are given by 

p(e , ) u 1 ft . n 

(V) 

r 0 o thorwi se 

where B. can be considered a ’’fitting bias” at E = Ej since 

E (6.) B. Var (€.) - 0. (8) 

Wliile all practical statistical properties of the regi-ession analysis become 
lost, the technique is not degraded as a numerical tool for apfjroxi mating with 
great precision complex fimctional forms over restricted regions with relatively 
simple functions. 
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To display the implicit linear form of <(>), we tahe the natural. Ip; 4 aritlim (it 
both sides of (6) mid obtain the equivalent torm 

In Y j In , I . I K. * , 

or 

y' , ; I • , E, , . ; (i 1. 2. ••• . tn) • (0) 

If nonnogativo weights . , , . . . , are available and the y, are eva.luatcd 

from y* = In h (K. ) = - I (cos .\,_j - cos E. ), then the well known result 



are the values of / ^ that minimize the weighted sum of squares 



1 1 


Tlius, b, = exp (bj) and are the classical v/eighted least squares estimates 
of/5,, ^2 respectively. 

The selection of the mesh points over an interval of Integration ! eJ, E‘; is 
obtained through an automated search approach utilizing the constrained 
weighted least squares process previously described. This technique selects 
subintervals of maximum size while retaining" a user selected error tolerance 
between the true and predicted function. 
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Furthermore, if one sets K* ~ U, K* “ , and derives ii sot of coefficients 
for the n selec ed subintervals in the manner previously defined, the functions 
h(K) and . (K) will have been fit for all values of E over which x has been as- 
sumed constant. This is true since h(K) is a periodic function with period 2 ^ 
and in the fundamental period (0, 2 ) is symmetrically distributed about the 
axis of symmetry F - ' . 

With the use of more general expression for density (c), the variational 
equations (1-17) integrated between E* and E* where 0 < E* < E* < , take on 


the following form: 
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E 

J-0 
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oxp 
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Having fit h(E) from (0 ), the integration of the variational equations between 

any two arbitrary limits E^ - 0, E* Ep takes the general form 


I ■— 1 ^ 

Aq^ : 2 ^ J ’ exp E} <1G^ (E) 

i - ^ * t 




exp [d'^_. E]} tlG^(E) 


’•E^tnofl 2 ^ 


rE^motl 2^ 

^ J 


exp (d". E) dG, (E) 


(13) 


t ’ t* exp (d^ E) dG^ (E) (k - 1. 2. •••. 6) 


Ey moti 2"^ 
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where rn *, ( ' = 1 , 2, . . 2n), are nonnegative integers; t , t', tl* (j -- 1, 2, - • 

n), (k = n H* 1, n + 2, . . 2n) are fui'.ctions of is an integer between 0 

and 2n - 2; is an integer betw'een 1 and 2n - 1; the grid points x ^ » ^,,,2 » 

. . defined from Xp , . x^; and (t*, t*d may be cither (t.r, t.), 

(tf, (tr, t! ), (tf, t ') depending upon the values K* mod 2 VI* mod 2 The 
logical structure of (13) and the determination of the above parameters is too 
lengthy to be given here but is presented in Appendix for completeness. Jt 
must be emphasized that (13) is general for all computations, but is valid only 
while X is assumed constant, whereby a new fit is obtained and new constanUi 
are determined to integrate (I- 13) by the general form (13). 

Using these results we can now solve the density' expression for the value 
of b within each of the selected or fitted intervals from the expression. 


b(L, U) r: 


~Pl ^ 

" ''l)^ 


(14) 


where the subscripts U and L still refer to the upper and lower points of an in- 
terval, For example, is that value of density' at the initial point of the inter-* 
val which Is known from the tables, and the corresponding density scale height 
can then be computed. The end point r^ is known as a function of E. Therefore is 
also ’advanced' in a similar fashion as during the fitting of the exponential. 

The above method of fitting King-IIele’s expression over several intervals of 
an orbit is sufficient to 'reproduce' the tabular model atmosphere to within 
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ninoly- eight percent at nil points, and sinmltaneotiTrty yield closed form intcj- 
grable equations for the variation of tfie dements. Table 1 provides a list of 
Cv')inputed versus static model atmospjiere densities for various altitudes and 
the differences, as experienced by the San Alarco-2 satellite during a revolution 
using the Spring-Tall Model with an Txospheric temperature of 1100"K. Where 
the density is high, as for example the critical region aroiind perigee, the differ- 
ence between calculated and tabular values are I-''”" grams/(km)^ or less. 

Since b(L, U) is determined by r^, to be at the upper end of each arc, then the 
density differences at this point become vanislringly small. 


III. SOLUI’ION OF VARIATIONAL I'QUATIONS 


In light of the analysis done in section II above, we now return to the vari- 
ational equations (1-1 S). If we now combine the radical terms in (1-18), we ob- 
tain a .^et of expressions containing the forms (1 ± x)"^, and (1 ± x)' ”, whex-e 
n = 1/2 and ^ = e‘ cos^ R. If one assumes that y will never get too large, 
that is, for drag satellites, e will not be larger than 0.2, then the above terms 
cai^ be expanded as, 


and 


. n(n - 1) o 

(1 i X)" ^ 1 ± nv + i 


O i X')'" - 1 V n X' i 


( 1 ) 


Adopting expi*cssion (II-l.l) for atmospheric density variation, the variational 
equations are then reduced to solving a set of indefinite integrals of the form 
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r*’’ . . , , . 

k I exp - {'x I{j ) Ej - <-'is H; ^ o s K s iir’ K ‘!K 

‘ 




where k, is simply tlie construU cociTicient of each corresponding iutegrui. It, 
on the other hand, tlie eccentricity is somewhat larger tlian normal, say around 
0.5 or so, then it will necessary to include several more terms in the ex- 
pansion (1) alxjve. This however is an algebraic problem, and compuUjUonally 
speaking has only the effect of changing the overall coefficients, { k ♦, in the 
final algorithm. The integrals are still of the form (2) above. Utilizing the fitting 
described in section IT, the exponential of the density can now be represented as 


bj(E) exp Ib^CE) 'El 




in which bj (P:) and are determined for one or several segments within an 

orbital revolution- Results indicate that b^(E) and b^iE) remain fixed for a 
considerable length of the satellite's lifetime. Refitting becomes necessary only 
when a and e change appreciably. This is found to occur more frequently near 
the end of the satellite's lifetime. In any event, the fitting procedure is 
instantaneous, and the calculation proceeds without interruption. 

Using these results^ equation (2) can now be rewritten as 

I bj(E) exp ib^CE) L] cos^ E*sin^ E dE (-1) 
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Here p takes the values 0, 1, and 2, Wlien p “ 0, ' takes on all values fron\ zero 
tlirough twelve. W^ien p 1, goes from zero through eleven, and wliere p - 2, 
ranges from zero through ten. V/lien p = 0, equation (4) is a recursion solu- 
tion giv'en b\* a sunmiation whose terms take the form 


exp [b (E)'E'' 


cos' E(b^ cos E + ^ s i n E) 

hfr^ 


4 . 


xp [b (E)*El cor/' "->EdE, 


where b^ is taken constant witJiin the subinteiwal (E^ , £^,1 of the interval [E*. E*l. 
When p - 1 or 2, the recursion solution summation terms become after some 
reduction, 

exp [b 2 (E)-E] cos^ EdE. (6) 

If in the drag acceleration one desires to express the velocity and density 
exponential so as to involve the product of the coefficient of the eii‘tn*s 
second zonal harmonic, and the air density / , one can follow the procedure of 
Sherrill (Reference 4) and write these terms as, 








and is approximated by = ^ 2 ^/ » where is the value of the earth’s 
equatorial I'adius (Reference 1). If is talvcn to be zero here, then tliese 

cross products in v are ignored and (7) and (6) reduce .o those cxpi'cssions 
given in Section I, Results from low perigee satellites such as OAR-901 (Cannon- 
ball) show that the c teiuns of equation (7) make differences only in the sixth 
decimal place in the variational equations. In addition, the cx’oss term in tlic 
exponential is 10^ Limes those corresponding terms in the velocity. Since this 
correction is accounted for rather simply by the fitting procedure, the cross 
term is retained in the exponential. 

The solutions of the variational equations are added to th^. epoch values of 
the Izsalv elements in the Vinti program, thus enabling the px'ogram to produce 
an ephemeris which includ es the effect of aerodynamic drag. In effect, the 
original separability is maintained, and the initial or boundary conditions are 
advanced using an analytic or closed form solution for the variational equations* 



IV. JU:SULTS 


In order lo lest the liMW theory, t\vo heav^^ air drag satellites with some- 
what dissimilar characteristics are considered. 'I'hese arc the Italian San 
>Larco-2 and the U.S. Air Force Caiuionhall (OAU-901). Data on these spacecrafts 
ifi as follows; 

San Marco- 2: 

Mass m = 129.27383 kgm. 

Projected Area A " 3425.3397 cm^. 

Drag Coefficient = 2.1 

Initial epoch; April 2G, 19G7, 10 hours, 12 minutes. 

Initial conditions: 

X = +0.58725272, x = -0.82890608 

y = +0.84923499, y = 0.56396878 

z = -0.05068537, z = +0.01219124 

Here, x, y, and z are in units of earth radii (6378.166 kilometers), and x, y, and 
z in units of earth radii per canonical iinit of time (806.812 seconds). 

Apogee height == 736.00 kilometers 
Perigee height *-= 205.60 kilometers 
Eccentricity = 0.0387 
Inclination = 2.87 degrees. 

Cannonball: 

Mass m =362.87392 kgm. 

Projected Area a = 3423.6195 cm^. 



Draj: Coefficient C,, = 2,1 

Initial Epoch: August 7, 1971, 0 hours, 20 ininutes. 

Initial Conditions: 

X = -0.9428339, x = -0.2G9775G5 

y = -0.28101G29, y = -0.40945775 

z ^ +0.28038380, z = -1.0092100 

Apogee height = 1957.20 kilometers 
Perigee height = 130.16 kilometers 
Eccentricity' = 0.1230 
Inclination = 92.00 degrees. 

In both cases, the boundary values of the atmospheric density profile 
given b> are taken from a static model, namely, the 196G U. S. Standard 
Air Force Supplements. The profiles used here include a Spring- Fa II 
model with an 1100 degree exospheric temperature, a Winter, 80U degree 
model, and a Summer 1000 degree model. While these profiles are 
adequate if chosen carefully, it is felt that a somewhat more sophisticated and 
dymamic model such as given by Jacchla (Uefcrence 5) would not only improve 
results but also make them more reliable. 

Figure 1, shows the variation over one oi'bital period (from time of in- 
sertion) of the semi-major axis, eccentricity, and (/• 2 + ^^^ 3 ) respectively for 
San Marco-2. Here, a^^ and are taken to be the initial values of these 
parameters. For the semi major axis, we have initially a secular decrease of 





108 parts in 10*' per i-ev'olution, with a |>oriodic variation superimposecl, of about 
25 parts in 10^' with the orbital peritxl. The eccentricity shows a secular de- 
crease of approxiniately 07 parts in 10^’ per revolution, and a pcricKJie variation 
superimposed on it with an amplitude of about 22 x 10"^, with the orl)ilal period. 
Figure 1 shows that ( , ‘ ^ ) shifts back :uid forth by about 8,3 seconds. Fi;;urc 2 

is a similar ^raph for Cannonball. It is seen that the variations for a and c 
here are somewhat larj^er than for San Marco-2 while the shift of ( - ^ is 

considci'ably less, I'his behavior is what one would expect considering the dif- 
ferences in the orbits. 

Figures 3 and 4 are graphs of semi- major axis versus time (days from in- 
sertion), for San Marco-2 and Cannonball during im actual lifetime study for the 
two satellites. l"ho circles are those values of a, predicted by the BMW theoiy 
using only the initial conditions given above, while the crosses indicate those 
arcs of data supplied by orbit improvement routines utilizing tracking or ob- 
servational datii to update the orbital elements. The orbital improvements were 
necessitated by the rapid tleterioration of the orbital parameter quality due to 
inaccurate force modeling (particularly air resistance accelerations) in the 
equations of motion. The accepted date of re-entry into the earth's atmosphere 
for San Marco-2 was on October 14, 19G7 at approximately 13 hours. Thus, the 
total lifetime was about 171 days* The program computed a lifetime of 

165 days, for a re-entry on October 8, 1967. Cannonball's re-entry date was 
approximately January 28, 1972, a lifetime of about 173 days. BMIV computed 




*\rr 


170 days. In diesc two cases, the pro^^rani evaluated the limits of the variational 
equations for values of the eccentric anomaly corresponding: to five day ijitervaJs. 
This was done so as to allow the BMW program to compare at intci'mediate 
points, its values for semi-major axis with the obseiwed ones. In actual prac- 
tice, these limits would be evaluated for those values of eccentric anomaly 
corresponding to those regions of an orbit over which a fit of the density variation 
remained fixed. Using a change in perigee height criterion (preselected at 1 km), 
results from Cannonball show that the first region is the first G5 days, the 
second is the next 55 days, luid so on, until within the Itist month of life the regions 
arc only of 5 day durations. As a result, the entire ephemeris is computed 
rapidly. Using the IBM 360/01 electronic digital computer*, the complete San 
Marco-2 and Cannonball ephenierides computed at 5 day intervals were exe- 
cuted in less than 18 seconds of computer CPU time. 

These results show that the liMlV progi*am was within 4 percent of the 
•’true’’ lifecimes of both satellites, despite the fact that Cannonball liad a re- 
latively high eccentricity and low perigee. In addition, it must be pointed out 
that the initial conditions used in the program were obtained from orbit im- 
provement routines other than the Vinti Orbit Determ* nation System (References 
6, 7), since this was the only data available. To be more consistent, and to Im- 
prove accrxacy, one should utilize, if possible, only those initial or epoch values 
obtained from the Vinti Orbit Determination System. 


In spite of the success vath San -Marco”2 amt Cminonball, ii7iprove- 

ments to the present BMW computer program are bein^. considered. I'rime 
among these is the incori>oration of a d^mamic mocfel atmospheric density pro- 
file such as Jaochia’s. As is well kno\vn a slight miscalculation in the exospheric 
temiTcraturc for a static model wouki drastically alter tlie resulting density 
profile, and consequently, the computed ephenioris of the satellite. 

The distinct advantage of the BMW method over numerical integration tech- 
niques is that being analytic or closed form, it is not subject to lai'ge eri'or 
accumulation due to roundoff and tnuKtation. In addition, an entire ephemeris 
can be obtained in a molter of a few seconds on present generation computers, 
while such a task might be prohibitive with a numerical integration. 

With tbo atmospheric approximation of this paper, equation (I1I-3) could be 
expressed in terms of Bessel functions, for example, 

exp (c cos E) cos” EdE .. I (c) 

- I 1 (c) 

= (c) -i I, (c) (n = 2) 

and so forth. The disadvantage here, besides the limited step size over which 
to perform the summation, is in the difficulty of handling intermediate points 



(n - 0) 
(n 1) 



within the limits. A Bessel function approach would appeiu' lo be better em- 
ployed in a study of atmospheric density inference from calculations of the 
period decide ment. 

King-Hole (Reference 3) has studied the contraction of orbits in a closed 
form manntir. Four situations are considered; 

1. Normal e. Phase 1: approx. 0.02 0.2 (3 » (ae/H^ ) ' 30 approx.) 

2. Normal e, Phase 2: 0 'e- 0.02 approx. (0 < (ae/lip)‘ 3) 

3. Circuhu’ orbits: e --- 0 (ae/Hp = 0). 

4. High eccentricity: e.: 0.2 (ae/Hp 30 approx.). 

It is felt that the regions covered by the equations of motion in these four cases 
are also covered by the BMW theory. Even though a static model atmosphere 
profile was used in the calculations, BMW does have the latitude of utilizing a 
density profile such as the Jacchia model. In addition, the reference orbit is 
not aji exact ellipse, but the Vinti orbit. As a result, it is felt that the BAfW 
method has contributed to the program for calculating both accurately and rapidly 
the orbits of satellites experiencing aerodynamic drag. 

At the present time, the Vinti differential correction algorithm (Reference 7) 
employing a classical weighted least squares technique is being appropriately 
modified to accommodate the augmented force model. The bulk of the imple- 
mentation requires the reformulation of the partial derivatives {^y. /^q^it^)} 
to reflect air drag where {y .} are tracking observables and {q. (t^)} are the 
set of ’’epoch** Vinti orbital elements and an atmospheric parameter iteratively 



estimated by the differentia. I correction technique. It appears that complete 
analyticity can be retained for these partial derivative expressions, thus rnerj^in^ 
the theoretical developments presented here with the practical application of 
orbit estimation. This work will ha forthcoming in a future report. 
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Figure 3. Variation of Semi- Major Axis for San Marco- 2 
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APPKNDIX 


DEVELOPMENT OF VAIUATIONAL EQUATION CONSTANT'S 

For k = 1, 2, • . 6; the general expression for the cimnge in the Vinti 

orbital elements due to atmospheric drag tetween the arbitrary limits E*, E* is 
given by; 


J ” It exp EJ } clGy (E) 

T2^-x 

'ZZ J * ‘'’‘P > dC^(E) 

i=0 


(A-1) 


^ Ej 5T1od 2*^ 

+ J exp [d* E] dCj^ (E) 


X 


t* exp [d* E3 dG, (E) 


E| mod 2“^ 


Without loss of generality, assume that the fitting process is perform over 
the interval [n,7Tj i.e. 

0 = X„ < Xj < Xj < ••• < X„_, < X„ =7T. 

the coefficients t^ are expressed by 





t. =c, LI +l)(L. U)x^ (cos - cos E)^l (( " 1. 2 - •• n) (A-2) 

where o is a previously defined fitted parameter in the interval 1- , x/i 
t = 1, 2, . . . n; X = ae, and b(L, U) is derived by (11-14). From syjnmetry 
considerations define 

’‘2n--e ^ (A -3) 

I in-1. J = <2;/ - x.^, 2t 7 _ X,f _ , ) (^-4) 

*^ 2 n-^+i ~ oxp(2n [l +b(L, U) x^ (cosEj - cos E)^] t - 1, 2. ... n (A-5) 

^ ^ (A-C) 


-- (A-7) 


where is the other previously defined, subinterval dependent, fitted parameter. 


The remainder of this appendix shows the procedure utilized to derive the 


values of the non negative integer multiplicity factors .U ^ = 1. 2, . . . n; 

the value ^ , an Integer behveen 0, 2n - 2 ; the value of , an Integer between 
1, 2n - 1; and (t^, tp which may be (t^, t^), (t^, t^), (t^, tp, or (t^, t^) dependlniT 
upon whether E* mod 2?/ and E* mod 2 t 7 are less than or equal to or greater 
than 77, for arbitrary E*, E* and for different modes of operation of the orbit 


£ 

if: 

jr-; 


A-2 



program. Prior to defining this procedure, consider the following definitions 
concerning (A) modes of operation of the Yinti orbit generator program, (B) 
integration intervals of the variational equatijns, and (C) Fortran variable 
definitions. 

(A) Mo<Je I 

The Vinti orbit generator program is operating in mode I when the time 
span of the variational equation integration [or implicitly, the integration 
interval (E*,E*)1 is sufficiently small so that there is not a complete fitting 
subinterval imbedded between the integration limits. Mode I will be the 
dominant mode for definitive ephemerides or an ephemeris computed to support 
a differential correction. 

Mode II 

The Vinti orbit generator is operating in mode n when the time span of 
variational integration has one or more imbedded fitting subintervals be- 
tween the limits of variational equation integration. This mode will be 
utilized in lifetime studies when large intervals of integration will be per- 
formed. Note that mode II implies the assumption of a valid fit (or x = 
ae constant) over large periods of time than does mode I where new fits will 
be performed as frequently as is required. 

(B) Class 1 Interval 

A class I interval is an intsrval such tliat both boundary points are adjacent 
grid points of the mesh x , x^. 



Class II Interval 


A class n interval is an interval such that one boundary point of the interval 
is a point x . of the mesh x^, x j, • . . , x„ while the other boundary point is a 
non mesh point between 

X . and X . , . 

3 3 + ' 


(C) Fortran Variable Definitions 
Fortran Variables 
El 
E2 

NOINTV 


Mathematical Description 


2n 




COEFF (I.J) J r. 1 


I=- 


MOINTV 


INOINTV + 1 


c , toe 


+ l ‘^2n 


I 2. NPiHY , 1 


J - 2 



d, tod„ 




2n 


J = 3 {1 = X, NOINTV + 1 Xq. X,. ••• 
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array (I, j) 


1=2 


1=3. NOINTV 


I = 3. 


j _ KOINTV 


J - 1 

E* mod 27] 

J = 2 


) 

0 or I 

J = 4 


J = 5 

V, 



I J = 2 Ej mod 1v 
J = 3 0 or 1 

J = 4 t* 

, J = 5 d* 


4 . 2 


{ 


J = 1 
J = 2 


’‘o. X,. 

X,. X,. 


nointv 



■n-t 


J = 3 

m n+>e 

J = 4 


J= s 

t 
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The following describes the subtasks associated with defining *.he full set of 
integration parameters in (A-1). Th^i detailed logic is preseni ^d in the flow- 
cliart given in Figure A-1. 

@ This algorithm will test whether E* E* are within the same sub- 
intervali will store the proper interval of integration and associated constants 
into ARRAY, and exit. This simplified logic Is most useful in Mode I orbit 
generations. 

@ This algorithm defines the Class n Intervals associated with the arbi- 
trary limits E*, E*, (i.e, [x^, E* mod 2 77 ] and [E* mod 2rr, x^j ) and stores them 
and the associated constants into ARRAY. 

(C) This algorithm checks to see if any of the Class II intervals are near 
77/2 or 3 7t/2 in order to avoid underflow difficulties during evaluation of the 
integrals. 

@ To assist the task of deleting unnecessary computations in Mode I gen- 
erations, this algorithm establishes the necessary intervals and associated con- 
stants when rf and E* are in adjacent subintervals and in the same multiple 
of 2r. 

@ This logic stores the entire set of Class I subintervals and a^^scciated 
constants into ARRAY and checks the Class n subintervals for negligible length. 
If such a Class 11 interval Is found, it is eliminated from consideration. -An 
arbitrary interval length of 1 x 10 ’^ is the criterion presently used. 


A-6 


^ - ■ 


This algorithm defines m^> m . It is predominantly uso<; in Mode 11 
orbit generations. 

@ This algorithm redefines ARRAY to eliminate Class! intervals with a 
multiplicity factor of zero and Class n intervals whose multiplicity factor has 
been reset from one to zero by when its length is negligible. 
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